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A time-dependent formulation for electron-hole excitations in extended finite systems, based on 
the Bethe-Salpeter equation (BSE), is developed using a stochastic wave function approach. The 
time-dependent formulation builds on the connection between time-dependent Hartree-Fock (TDHF) 
theory and configuration-interaction with single substitution (CIS) method. This results in a time- 
dependent Schrodinger-like equation for the quasiparticle orbital dynamics based on an effective 
Hamiltonian containing direct Hartree and screened exchange terms, where screening is described 
within the Random Phase Approximation (RPA). To solve for the optical absorption spectrum, 
we develop a stochastic formulation in which the quasiparticle orbitals are replaced by stochastic 
orbitals to evaluate the direct and exchange terms in the Hamiltonian as well as the RPA screening. 
This leads to an overall quadratic scaling, a significant improvement over the equivalent symplectic 
eigenvalue representation of the BSE. Application of the time-dependent stochastic BSE (TDsBSE) 
approach to silicon and CdSe nanocrystals up to size of « 3000 electrons is presented and discussed. 


I. INTRODUCTION 

Understanding electron-hole excitations in large molec¬ 
ular systems and nanostructures is essential for develop¬ 
ing novel optical and electronic devices.^® This is due, 
for example, to the exponential sensitivity of the photo¬ 
current characteristics to the excitonic energy levels and 
the sensitivity of the device performance to the optical 
oscillator strength. It becomes, therefore, a necessity to 
develop accurate theoretical tools to describe the exci¬ 
tonic level alignment and the absorption spectrum, with 
computational complexity that is scalable to systems of 
experimental relevance (thousands of atoms and more). 

There is no doubt that time-dependent density func¬ 
tional theory (TDDFT)® has revolutionized the field of 
electronic spectroscopy of small molecular entities.^^®^ 
TDDFT provides access to excited state energies, ge¬ 
ometries, and other properties of small molecules with a 
relatively moderate computational cost, similar to config¬ 
uration interaction with single substitutions (CIS) in the 
linear response frequency-domain formulatioii^ (O , 
where N is the number of electro^, or even better us¬ 
ing a real-time implementatioEpSlMl (jj In prin¬ 

ciple TDDFT is exact but in practice approximations 
have to be introduced. The most common is the so-called 
time-dependent Kohn-Sham (TDKS) method within the 
adiabatic approximation, which has been applied to nu¬ 
merous challenging problem^i^MS] .^^^h great success. 
However, TDKS often fails,particularly for charge- 
transfer excited states, multiple excitations, and avoided 
crossings. In the present context, perhaps the most sig¬ 
nificant failure of TDKS is in the description of low-lying 
excitonic states in bulk.l22HIZI 


An alternative to TDDFT, which has mainly been ap¬ 
plied to condensed periodic structures, is based on many- 
body perturbation theory (MBPT). The most com¬ 
mon flavors are the GW approximatiorP^ to describe 
quasiparticle excitations (G indicates the single-particle 
Green function and W the screened Goulomb interac¬ 
tion) and the Bethe-Salpeter equation (BSE))^ to de¬ 
scribe electron-hole excitations. Both appr oaches of¬ 
fer a reliable solution to quasiparticle^^HlIl and opti- 

Ca l53 | 54 | 56 | 57 | 661 79lMl excitations, even for situations where 
TDKS often fails, for example in periodic systemd^^HIZIZll 
or for charge-transfer excitations in molecules.^ How¬ 
ever, the computational cost of the MBPT methods is 
considerably more demanding than for TDKS, because 
conventional techniques require the explicit calculation of 
a large number of occupied and virtual electronic states 
and the evaluation of a large number of screened ex¬ 
change integrals between valence and conduction states. 
This leads to a typical scaling of O (A^®) and limits the 
practical applications of the BSE to small molecules or 
to periodic systems with small unit cells. 

Significant progress has been mad e by c ombining ideas 
proposed in the context of TDDFfrfsnilT] and techniques 
used to represent the dielectric functiorP^ based on den¬ 
sity functional perturbation theory.!^ This leads to an 
approach that explicitly requires only the occupied or¬ 
bitals (and not the virtual states) and thus scales as 
O X X where Nk the number of points in 

the Brillouin zone and Ng the size of the basis. Even with 
this more moderate scaling, performing a Bethe-Salpeter 
(BS) calculation for large systems with several thousands 
of electrons is still prohibitive. 

Recently, we have proposed an alternative formula- 
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tion for a class of electronic structure met hods ranging 
from the density functional theory (DFT)P^*25| M0ller- 
Plesset second order perturbation theory (MP2)f2^*2Zlthe 
random phase approximation (RPA) to the correlation 
energy}^ and even for multiexciton generation (MEG) 
But perhaps the most impressive formulations are that 
for calculating the quasiparticle energy wit hin the GW 
many-body perturba tion correction to DFTR22I and for 
a stochastic TDDFTP^ The basic idea behind our for¬ 
mulation is that the occupied and virtual orbitals of the 
Kohn-Sham (KS) Hamiltonian are replaced by stochas¬ 
tic orbitals and the density and observables of interest 
are determined from an average of stochastic replicas in 
a trace formula. This facilitates “self-averaging” which 
leads to the first ever report of sublinear scaling DFT 
electronic structure method (for the total energy per elec¬ 
tron) and nearly linear scaling GW approach, breaking 
the theoretical scaling limit for GW as well as circum¬ 
venting the need for energy cutoff approximations. 

In this paper we develop an efficient approach for cal¬ 
culating electron-hole excitations (rather than charge ex¬ 
citations) based on the BSE, making it a practical and 
accessible computational tool for very large molecules 
and nanostructures. The BSE is often formulated in the 
frequency domain and thus requires the calculation of 
screened exchange integrals between occupied and virtual 
states. Instead, we introduce concepts based on stochas¬ 
tic orbitals and reformulate the BSE in the time-domain 
as means of reducing GPU time and memory. The real¬ 
time formulation of the BSE delivers the response func¬ 
tion (and thus the optical excitation spectrum) with¬ 
out requiring full resolution of the excitation energies, 
thereby reducing dramatically the computational cost. 
This is demonstrated for well-studied systems of sili¬ 
con and GdSe nanocrystals, covering the size range of 
N « 100 — 3000 electrons. Within this range, we show 
that the approach scales quadratically (O {N^)) with sys¬ 
tem size. 


II. THEORY 

In this section we review the symplectic eigenvalue for¬ 
mulation of the BSE and then build on the connections 
between configuration interaction with single substitu¬ 
tion (GIS) and time-dependent Hartree-Fock (TDHF) to 
formulate a time-dependent wave-equation for the BSE. 


A. Symplectic Eigenvalue Bethe-Salpeter Equation 

Within linear response, one can show that the BSE 
is equivale nt to solving the symplectic eigenvalue prob- 
lenpmni 



where 

^ ^ (-4 (^) 

with 

A = D + 2K^ + 

B = 2K^+K^. (3) 

The diagonal {D), exchange {K^) and direct {K'^) terms 
are given by (we use f,j, and k... as occupied (hole) 
state indices, a, b, and c... as unoccupied (electron) 
states indices, and r, s, and t... for general indices): 


Dia^bj 

— (^a ^i) babbij 


(4) 

-^ia,bj 

= {4’a4>i\Pc\(l>b4’j) = 




(r) (pa (r) vc{\y 

(r')</>b (r') 

(5) 

-^^ia^bj 

= {(l>a(t)b\W\4>i(t)j) = 

If Mr’ 



(r) (j)a {r)W (r. 

{y')4>, (r'). 

(6) 


Here, Sa and Si are the quasi-particle energies for the vir¬ 
tual and occupied space (which can be obtained from a 
DFT+GW calculation or from an alternative suitable ap¬ 
proach) and (j)a (r) and (j)i (r) are the corresponding quasi¬ 
particle orbitals; vc is the Goulomb potential while W 
is the screened Goulomb potential, typically calculated 
within the Random Phase Approximation (RPA), which 
can be written in real space as: 


W (r, r', 0) = nc (|r - r'|) + (r, r', 0), (7) 

with 


0)= JJ dY"dv'"vc (|r - r"|) 

^RPA (r",r"',0) (nc (|r'" - r'|) + fxc[Y'")5{v'" - r')). 

( 8 ) 

Here, fxci^) is the DFT exchange-correlation potential 
(if DFT is used to obtain the RPA screening, otherwise 
set /xc(r) = 0), and (r, r', 0) is the half-Fourier 

transform (at w = 0) of the real-time density-density 
correlation function within the RPA (the latter can be 
also obtained from TDDFT, as further discussed below). 
We note in passing that often the above is solved within 
the Tamm-Dancoff approximation (TDA),i ^^P°^l which 
sets 5 = 0 and thus only requires the diagonalization of 
the matrix A. 
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B. Time-Dependent Bethe-Salpeter Equation 
(TDBSE) 

The time-dependent formulation of the BSE fol¬ 
lows f rom the con nections made between CIS and 
TDHF.^S2110 3 |io7 | io 8| solving the TDHF equa¬ 
tions = fiHF (t) <l)j {v,t) for the occupied 

orbitals is identical to solving the symplectic eigen¬ 
value problem of Eq. 0 with (IIF(r,r',0) = 0. 

Here, hnF = i + Vion + vh {t) + kx (t) is the 
Hartree-Fock (HF) Hamiltonian, t is the kinetic en¬ 
ergy, Vion is the external potential, vh'>P (r) = 

/ dr'uc (|r — r'l) u (r', t) i/) (r) is the Hartree potential, 
and kx{t) tp{r) = -i / dr'p (r, r', t) uc (|r - r'|) V' (r') 
is the non-local exchange potential. n (r, t) = 
and p(r,r',t) = 

are the time-dependent electron density and density ma¬ 
trix, respectively. The connection to CIS is made by 
realizing that for 5W (r,r',0) = 0 and setting B = 0 
(the TDA), the symplectic eigenvalue problem of Eq. 0 
is nothing else but the CIS Hamiltonian. Thus, TDHF 
within the TDA and CIS are identical. 

We follow a similar logic and derive an adiabatic time- 
dependent BSE: 

d(j)^ (r, t) 

0 ( 9 ) 

where 7 is a perturbation strength (i.e., 7 = 0 is the 
unperturbed case, see Eq. (HI) with a screened effective 
Hamiltonian given by: 

hls = Kp + Vh (t) - Vh (t) + k]x i*) “ k°x {*) ■ (10) 

Here, hqp is the quasi-particle Hamiltonian which is 
typically determined from a GW calculation correcting 
the quasiparticle energies and orbitals of the underly¬ 
ing DFT. The GW approximation to hqp is rather diffi¬ 
cult to implement since it involves a non-local, energy- 
dependent operator. An alternative is to use a DFT ap¬ 
proach that pro vides an accurate description of quasipar¬ 
ticle excitations .(l^ninni However, since the exact model 
for hqp is not the central target of the present work, we 
represent it by a sim ple semi-empirical local Hamiltonian 
of the formffimilll 


i-qp 


' t + u 


psi 


( 11 ) 


where, as before t is the kinetic energy and Vps = Va 
is the empirical pseudopotential, given as a sum of atomic 
pseudopotentials which were generated to reproduce the 
bulk band structure, providing accurate quasi-particle 
excitations in the bulk. The semiempirical approach has 
been successfully applied to calculate the quasi-particle 
spectrum of semiconducting nanocrystals of various sizes 
and shapes.EanMIliniHI^ 

Eq. v]j{t)'ip{r) = 


In 


/ (|r — r'l) (r^, t) ^ (r) is the Hartree potential 

with n'^(r,f) = 2 |<?^J (r, and k]^ (t) tp (r) = 

— dr'p'’’(r,r',t) (r,r',0) Ip (r') is the screened 

exchange potential with (r, r', 0) given by Eqs. Q 

and ^ and p^(r, r',t) = (r,i)- The 

application of kj-^ip (r) is further discussed below. 



Figure 1. Comparison of BS calculations using the symplectic 
eigenvalue (Eq. 0 . a:-symbols) with the frequency dependent 
dipole-d ipol e correlation generated from TDBSE (Eqs. (10 1 , 
H and ( |11[ ), solid lines) for SiH 4 . Black: TDBSE with 7 = 0 
(TDHO) compared with eigenvalues of Eq. 0 setting and 
RV’ to zero. Red: TDBSE with fpgg = hqp + vjj (t) — Vh (t) 
(TDH) compared with Eq. ||^ setting K'^ to zero (TDH). 
Green: TDBSE with 7 = 10“ au (TDBSE) compared with 
Eq. ||^, both for e = 5. 


In analogy with the relations derived between TDHF 
and its eigenvalue representation, it is clear that the time- 
dependent formulation for the BSE given by Eqs. 0 and 
(101 is identical to the full symplectic eigenvalue prob¬ 
lem of Eq. 0. In Fig. [T] we compare the results for 
SiH 4 on a 8 X 8 X 8 grid generated by propagating the 
occupied orbitals with the Bethe-Salpeter Hamiltonian 
(10 1 (TDBSE) to the exact diagonalization of Eq. 0 
(static approach). We use a local semi-empirical pseu¬ 
dopotential that has been applied successfull y to study 

the optical properties of silicon nanocrystals 

For both the direct approach and the TDBSE we ap¬ 
proximate W(r,r',0) by (|r — r'|), where e is a 

constant screening parameter. The idea is to confirm 
that the eigenvalues of Eq. 0 and the time-dependent 
version of the BSE are identical (validating both the the¬ 
ory and the implementation). 

The time-domain calculations are based on a linear- 
response approach to generate the dipole-dipole corre¬ 
lation function d (t) and its Fourier transform d (w) = 
dte^‘^*d pt). In short, we perturb the occupied eigen¬ 
states {(pj (r)) of hqp at f = 0: 


(p] (r,t = 0) = e (r), 


( 12 ) 
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where for simplicity, we assume that the dipole is in the 
^-direction. We then propagate these orbitals accord¬ 
ing to Eq. and generate the dipole-dipole correlation 
function: 


d{t) = — J dr z {v? (r, t) — tP (r, t)) , (13) 

where as before n'>'(r,t) = 2 ^^ \(f- (r,t)|^and 7 is a 
small parameter representing the strength of the pertur¬ 
bation, typically 10 “^ — . 


The agreement for the position of the excitations (solid 
lines) generated by the time-domain BSE is perfect with 
the static calculation (x-symbols), as seen in Fig.[^ The 
resolved individual transitions are broadened reflecting 
the finite propagation time used for the time-domain 
calculations. We find that in some cases the oscillator 
strength is very small and thus a transition is not ob¬ 
served in d (w). 



Figure 2. Average momentum along the z-direction calculated 
in two ways (see text for more details) for SiH 4 using the 
TDH (black curves) and the TDBSE (red curves) methods. 
Solid and dashed curves where generated using the expecta¬ 
tion value of the momentum (Eq. @) and the numerical 
derivative of the expectation value of the position (Eq. (dt). 
respectively. Inset: same for longer times. 


III. TIME-DEPENDENT STOCHASTIC 
BETHE-SALPETER EQUATION 


An additional important test of the TDBSE formal¬ 
ism is whether the Hamiltonian in Eq. (10 1 preserves the 


Ehrenfest theorem (see Appendix B for more details). 
Naturally, this would be the case if hqp would include 
the terms (t) and (t), such that they cancel out 
for h'^g. However, for an arbitrary choice of hqp this needs 
to be confirmed. In Figure|^we plot the average momen¬ 
tum for SiH 4 calculated in two different ways. The solid 
curves were obtained directly from: 


(P W) 




dr(j)] {r,t)* ^4>] (r,t), 


(14) 


while the dashed curves were obtained by taking the nu¬ 
merical time derivative (central difference) of the expec¬ 
tation value of r it) : 

^ (q (0)=2^ E / y • 

( 15 ) 

The agreement is not perfect but improves with decreas¬ 
ing the time step 5t (not shown here). We also show 
the results for the time-dependent Hartree (TDH), i.e., 
ignoring the screened exchange term in hP^g. The devia¬ 
tions observed for TDBSE and TDH are similar, although 
for TDH the Ehrenfest theorem holds exactly and thus 
the agreement should be perfect. The difference are as¬ 
sociated with numerical inaccuracies resulting from the 
finite time step and grid used in the calculation. The 
inset shows that the deviations are insignificant even at 
much longer times over many periods. 


We consider two formulations for the time-dependent 
stochastic BSE (TDsBSE). The first approach is a direct 
generalization of the appro ach we have recently devel¬ 
oped for the stochastic TDH,Ii 2 l] in which we describe an 
efficient way to account for the the screened exchange 
term in the h'^g. This approach works well for short 
times, however, unlike in TDH, the inclusion of an ex¬ 
change term requires an increasing number of stochastic 
orbitals with the system size. The second approach of¬ 
fers access to timescales relevant for most spectroscopic 
applications at a practical quadratic computational cost. 


A. Extending the Stochastic TDH to Include a 
Screened Exchange Term 

We limit the discussion, in the body of this paper, to 
the case where W (r, r', 0) is replaced by e~^vc (|r — r^l), 
where e is a function of |r — r'|. The algorithm for the 
TDsBSE is based on the following steps: 

1. Generate Nq stochastic orbitals Q (r) = 

e®G(r) j^SV, where (r) is a uniform ran¬ 
dom variable in the range [ 0 , 27r] at each grid point 
(total of Ng grid points), SV is the volume element 
of the grid, and j = 1,..., N^. The stochastic 
orbitals obey the relation 1 = (|C)(CI)^ where 
(• • • )^ denotes a statistical average over (). 

2. Project each stochastic orbital Q (r) onto the oc¬ 
cupied space: \^j) = \ dfi [l^ — hqp] \C,j) , where 
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Op {x) = ^erfc {f3 (x)) is a smooth representation of 
the Heaviside step functioiPH and /r is the chemical 

potential. The action of is performed using a 
suitable expansion in terms of Chebyshev polyno¬ 
mials in the static quasi-particle Ha milt onian with 
coefficients that depend on and /3.E^ 

3. Define non-perturbed and perturbed orbitals for 

t = 0 to the orbitals: ^9(r,t = 0) = (r), 

(r,t = 0) = (r). For the absorption 

spectrum, the perturbation is given by v (r) = 
and a = x,y, z. 

4. Propagate the perturbed (^J(r,t)) and unper¬ 
turbed (r, t)) orbitals according to the adiabatic 
time-dependent BSE: 

Qfl (v t') ^ 

^ (i) Cj (r, t). (16) 

Use the split operator technique to perform the 
time propagation from time t to time t + At: 

ft! g“27i 

X g“ 2xf^*g“ K {^Ix (*))‘^* 

where propagator step involving the non-local 
screened exchange is applied using a Taylor series 
(in all applications below we stop at ): 

H (^A ft! 

1-^ (0 - fc°x (o) At -h • • • (18) 

5. The application of h'^g (t) is done as follows: 

(a) The kinetic energy is applied using a Fast 
Fourier Transform (FFT). 

(b) The Hartree term is generated using convolu¬ 
tion and FFT with the density obtained from 
the stochastic orbitals: 

■ (19) 

(c) The time-consuming part of the application 

of h'^g on a vector ip in Hilbert space is 
klx (t) — (0- This operation scales as 

O {NNgrid) and one needs to carry this for all 
occupied states, leading a O (N^Ngrid) com¬ 
putational scaling. To reduce this high scaling 
resulting from the exchange operation we use 
the same philosophy underlying this work, i.e., 


replacing summation with stochastic averag¬ 
ing. In practice we therefore replace the sum¬ 
mation over occupied orbitals in the exchange 
operation by acting with very few rig <?; 
typically rig = 1 — 16, stochastic orbitals write 
the exchange operation as: 

Kx (^) V' (r, i) = — 51 (’’’ ^) 

X=1 

xfdr'e~^vc (|r - r'|) 77^ (r',t)* ip (r', t). ( 20 ) 

The key is that these stochastic orbitals are 
defined as a different random combination of 
the full set of orbitals at any given time stepr;^ 
are defined as random superpositions of the 
stochastic orbitals: 

1 

= (21) 

To improve the representation of the stochas¬ 
tic exchange operators, the random phases 
axj (t) are re-sampled at each time step. Note 
that the same phases are used for both rjj (r, t) 
and ?7°(r,f). This use of stochastic orbitals 
reduces the overall scaling of the method to 
quadratic, since rig does not dependent on the 
system size. 

In Fig. we show the calculated d (t) and S (t) = 
/p dsd{s) for a series of silicon nanocrystals. We used 
riQ = 16 which leads to results that are indistinguishable 
from = A^ (though even a smaller would have been 
sufficient). We used a constant value fore = 5 and the 
time step was At = 0.025au. 

In general, we find that the results converge up to a 
time Tc and then the signal diverges exponentially. Sev¬ 
eral conclusions can be drawn from these calculations: 

1. The stochastic approximation to d{t) oscillates 
about zero up to a time tq, but this is followed 
by a gradual increase which eventually leads to di¬ 
vergence (upper panels of Fig. [^. 

2 . Tc increases with the number of stochastic orbitals, 
A^, roughly as tc oc A“ with a = 1 — 2 (right 

panels of Fig. 1^. This is somewhat better than the 

1/2 

case for TDH for which tq roughly scaled as A^ . 

3. Tc decreases with increasing system size roughly as 
^, where Ae is the number of electrons (left panels 
of Fig. 1^. Therefore, to converge the results to a 
fixed Tc one has to increase N(^ roughly linearly 
with the system size . This leads to a quadratic 
scaling of the approach. In TDH the opposite is 
true, Tc increas es w ith increasing system size due 

to self-averaging.IiSI] 
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time (au) time (au) 


Figure 3. Upper left: Dipole-dipole correlation function (Eq. normalized to the number of silicon atoms in the nanocrystal 
(Nsi) for several nanocrystals sizes. For each size we use a different number of stochastic orbitals. Lower left: Integrate dipole- 
dipole correlation S (t) = dsd (s)^. The onset of divergence scales roughly linearly with the size. Upper right: Dipole-dipole 
correlation function normalized to the number of silicon atoms for SisrHye for different values of N^. Lower right: Corresponding 
values for S (t). 


4. To reach times sufficient for most spectroscopic ap¬ 
plications, the number of stochastic orbitals ex¬ 
ceeds that of occupied states {N(^ > Nocc)- 

To conclude this subsection, we find that this version of 
a TDsBSE scales roughly quadratically with the system 
size, rather than sub-linearly for TDH. Furthermore, to 
calculate the response to meaningful times, the naive ex- 
tention of the TDH to include exchange requires a rather 
large number of stochastic orbitals (fV,^), often much 
larger than the number of occupied orbitals. However, 
it is sufficient to represent the operation of the exchange 
Hamiltonian with a relatively small set of linear combi¬ 
nation of all stochastic orbitals (n^;). We next show how 
the method can be improved significantly increasing tq 
to values much larger than required to obtain the spec¬ 
trum in large systems. 


B. Time-Dependent Stochastic Bethe-Salpeter 
with Orthogonalization 

To circumvent the pathological behavior observed 
above, we propose to orthogonalize the projected stochas¬ 
tic orbitals (after step “2”). This requires that be 
equal to the number of occupied states Nocc- However, 
this makes the TDsBSE stable for time-scale exceeding 
50fs, which for any practical spectroscopic application for 
large systems is more than sufficient. Formally, since the 
number of stochastic orbitals (equal to the number of 


occupied states) increases linearly with the system size, 
the approach scales as O {N/^Ng). The orthogonalization 

step scales formally as O (^N^Ng'^, however, for the size 
of systems studied here, it is computationally negligible 
compared with the projection and propagation steps. 

In Fig. we compare the dipole-dipole correlation 
function computed from the TDsBSE with = 1 to 
the direct TDBSE approach for silicon nanocrystals of 
varying sizes (SiasHae, SisrHyg, Sii 47 Hioo, SiaaaHigg, 
and SiyoaHaoo)- The purpose is to demonstrates the 
power of the TDsBSE approach with orthogonaliza¬ 
tion. Therefore, for simplicity W (r, r', 0 ) is replaced by 
(|r ~ r'l) with e = 5 for all system sizes. Clearly, 
even when = 1, the TDsBSE is in perfect agreement 
with the direct TDBSE approach. The cubic scaling of 
the later limits the application to small NCs or to short 
times. 

In Fig. 1^ we plot the TDsBSE absorption cross sec¬ 
tion (cr (w) = f drdr' zx (r, r',a;) z') compared to 
the absorption cross section computed by ignoring the 
electron-hole interactions for a wide range of energies. It 
is practically impossible to obtain the absorption cross 
section over this wide energy range by a direct diag- 
onalization of the symplectic eigenvalue equation (cf., 
Eq. 0). Thus, so far the BSE has been applied to rel¬ 
atively small nanocrystals or by converging only the low 
lying excitonic transitions, even within the crude approx¬ 
imation where W (r, r', 0 ) is replaced by e~^vc (|r — r^D- 
As far as we know the results shown in Fig. are the 
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Figure 4. The dipole-dipole correlation function calculated using the TDsBSE approach with orthogonalization (red curves) 
compared with a direct time-dependent BSE approach (black curves). Note that the direct (i.e., non-stochastic) BSE approach 
is so expensive due to the full-exchange operation that it was not done for the largest NCs and was only followed for short 
times for intermediate size NCs. 


first to report a converged BS calculations for NCs of 
experimentally relevant sizes. We used a constant e in 
each run, with values of 5, 6.2, 7, 8.2, and 8.8 taken from 
Ref. ll2d| for the silicon NCs (in ascending order) and 4.5, 
5, 5.2 and 5.4 for the CdSe NCs taken from Ref. 11131 The 
inclusion of a more accurate description of the screening 
as proposed in detailed in Appendix A is left open for 
future study. 

For both types of NCs there is a shift of the on¬ 
set of absorption to lower energies with increasing NC 
size due to the quantum confinement effect. The ab¬ 
sorption cross section of the smallest NCs is character¬ 
ized by detailed features, which are broadened and even¬ 
tually washed out as the NC size increases. For sili¬ 
con NCs, the semi-empirical pseudopotential model over¬ 
emphasizes the lowest excitonic transition in comparison 
to the p lasmonic r esonance observed at ~ lOeV using 
TDDFT.I^^lSlMIIoIl jjijgggg tjie gp^j^ Qf lowest 

excitonic peak observed experimentally for bulk silicon 
and reproduced by the BSE approach,l2^1t^ but not by 
the current model ignoring electron-hole correlations .1^^ 
The fact that the current calculation does not capture 
this split could be a consequence of the approximation 
used to model the screening. 

The results for silicon NCs seem to imply that the in¬ 
clusion of electron-hole interactions leads to a blue shift 
in the absorption cross section (black curve is shifted to 
higher energies compared to the red curve). Since silicon 
is an indirect band gap material, the onset of absorption 
is not a good measure of the strength of the electron- 
hole interactions. Indeed, when the approach is applied 
to CdSe NCs (lower panels of Fig. the inclusion of 
electron-hole interaction clearly shifts the onset of ab¬ 
sorption to lower energies. 


IV. CONCLUSIONS 

We have developed a real-time stochastic approach to 
describe electron-hole excitations in extended finite sys¬ 
tems based on the BSE. Following the logic connect¬ 
ing TDHF and CIS, we showed that a solution to a 
Schrodinger-like time-dependent equation for the quasi¬ 
particle orbitals with an effective Hamiltonian contain¬ 
ing both direct and screened exchange terms is equiv¬ 
alent to the symplectic eigenvalue representation of the 
BSE. A direct solution of the TDBSE leads to at least 
cubic scaling with the system size due to the need to 
compute all occupied quasiparticle orbitals and the com¬ 
plexity of applying the screened exchange term to pre¬ 
form the time propagation. The lower bound is similar 
to the scaling of the TDHF method and thus, limits the 
application of the TDBSE approach to relatively small 
systems. To overcome this bottleneck, we developed a 
stochastic approach inspired by our previous w ork on 
stochastic (sGW) and stochastic TDDFT,ff2il in 

which the occupied quasiparticle orbitals were replaced 
with stochastic orbitals. The latter were then used to 
obtain both the RPA screening using the approach de¬ 
veloped for the screening in sGW and the exchange po¬ 
tential by extending the approach used todescribe the 
Hartree term in TDsDFT. Both the RPA screening the 
application of the exchange potential scale nearly lin¬ 
early with system size (as opposed to quadratic scaling 
for example for the exchange potential). The number of 
stochastic orbitals required to converge the calculation 
scales with system size and thus, the overall scaling of 
the TDsBSE approach is quadratic (excluding the cubic 
contribution from the orthogonalization of the stochas¬ 
tic orbitals, which for the system sizes studied here is a 
negligible step). 

We have applied the TDsBSE approach to study opti¬ 
cal excitations in a wide range of energies (up to 30 eV) 
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Figure 5. Upper panels: The absorption cross section for silicon NCs computed by Fourier transforming the TDsBSE dipole- 
dipole correlation function (black curves) and the corresponding absorption cross section computed for hqp, i.e., by ignoring 
electron-hole interactions. 


in silicon and CdSe nanocrystals with sizes up to « 3000 
electrons (« 3 nm diameter) and compared the re¬ 
sults with the quasiparticle excitation spectrum obtained 
within the semi-empirical pseudopotential approach. For 
both systems, we find that including electron-hole cor¬ 
relations broadens the spectral features and shifts the 
oscillator strength to higher energies due to amplifica¬ 
tion of a plasmon resonance near 10 eV. For silicon we 
find a surprising result where the onset optical excita¬ 
tions seem to shift to higher energies compared to the 
quasiparticle excitations. This is a result of two factors. 
First, silicon is an indirect band gap material and th the 
onset of optically allowed transitions is above the the 
lowest excitonic state. Second, the inclusion of electron- 
hole interactions via the BSE leads to an amplification 
of a plasmon resonance at « 10 eV shifting the oscillator 
strength to higher energies at the expense of the lower 
frequency absorption. These combined effects lead to an 
apparent shift of the absorption onset to higher energies 
when electron-hole interactions are included. This is not 
the case for CdSe, where the onset of optical excitation is 
below the onset of the quasiparticle excitation, as expect 
for a direct band-gap material. 

The TDsBSE provides a platform to obtain optical 
excitations in extended systems covering a wide energy 


range. To overcome the divergent behavior at long times, 
it is necessary to increase the number of stochastic or¬ 
bitals as the size of the system increases. We are working 
in improvements of this flaw and if solved, an even faster, 
linear scaling BS approach will emerge. This and other 
improvements as well as more general applications will 
be presented in a future work. 
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Appendix A: RPA screened exchange for TDsBSE 

The above approach assumes that lF(r, r',0) = 

uc(|r —r'l) + 5Wrpa{y,v' ,0) is approximated by 
(|r — r'l). In typical BS applications, one uses the 
RPA screening to describe W (r, r', 0) = uc (|r — r^|) + 
(51F (r, r', 0). The stochastic formalism, however, fur¬ 
nishes a potentially viable approach to overcome the as- 
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sumption made to obtain W (r, r', 0) in this work. In the 
linear response limit, SWupa can be written as: 

5Wb.pa (r, r', ^) = JJ dr"dr'"vc (|r - r"|) 

XRPA (r",r"',0) X 

(uc(|r"'-r'|) + /xc(r'W"-r')), 

( 22 ) 

and we are concerned with the application of k^x (t) on 
Ip or more accurately, the portion that depends on 
the screening: 


tic time-dependent equations: 


('^) X] (r, t) , (27) 

Here, one can take h]^px (t) = hqp or hppA (t) = 


For the latter case, vhxc W (r) = / dr' + 
vxc (i")) E'-iid vxc (n (r)) is the local density (or 
semi-local) approximation for the exchange corre¬ 
lation potential. The density is obtained as an av¬ 
erage over the RPA stochastic orbital densities: 


i*) = (r, t) JdrSWppA (r, r', 0) r]2 (r', t)* ip (r', t). 


We first insert Eq. ® into Eq. 


(23) 


Nrpa 


^RPA 




i=i 


(28) 


5kPx {t) Ip (r, t) = r]2 (r, t) jj /dr'dr"dr'" 

VC (|r - v"\)xRPA (r",r'",0) {vc (|r'" - r'|) 

+ fxc{r"')5{v"' -v'))nZ{vPtr iPivPt). 

(24) 


Define a perturbation potential 

v~^ (r,t)= J dr' (uc (|r"'- r'l) 

+ fxc{r'")S{r'" - r')) r,2 (r', t)* ip (r', t) (25) 


4. Generate AnppA (r, r) = 

^ {v-RPA (’'>''”) ~ (r, r)^ and its half Fourier 

transformed quantity AhppA (r, 0 ) at w = 0 . 

5. Obtain the action of Sk^x i^) = 

rj2 (r, t) ff dr'dr"vc (|r - r'|) xrpa (r', r", 0) v'r (r", t) 
from SkfX (t) ip (r, t) = 

V2 (r, t) Jf dr'dr"VC (|r - r'|) AUrpa (r', 0). 

Step 1-5 need to be repeated at each time step At of the 
TDsBSE propagation. 


and rewrite Eq. (241 as: 


Appendix B: Ehrenfest theorem 


SkZx (t) Ip (r, t) = r]2 (r, t) JJ dr'dr" 

VC (|r - r'DxRPA (r', r", 0) v'^ (r", t). (26) 

The action of xrpa (r', r", 0) on v'^ (r", t ) is manageable 
by using a stochastic TDDFT algorithm)!^ 


Ehrenfest theorem asserts that a correct propagation 
must preserve the relation 


(q(t)) = i fiBsA ) 
For a TDBSE this relation is given by 


(29) 


1 . TakeNppA projected stochastic orbitals from the 
generated above. If Nrpa > N(^ generate ad¬ 
ditional projected stochastic orbitals following the 
prescription given in 1 and 2 above. This needs 
to be done just once, i.e., at the beginning of the 
calculation, generate enough projected stochastic 
orbitals to be used throughout the calculation. 


fiBsA 


(P (t)) 


m 


klx (t) - kex (t) A ) (30) 


where kJx {t) '‘P (^) = ~^ f dr'p"' (r, r', t) (r, r', 0 ) ip (r'). 

klx (t) - Kx (t) > q 


To satisfy the Ehrenfest theorem 
should vanish. The commutator of the exchange operator 
is given by: 


2. Apply a perturbation at r = 0: x] = 0) = 
(r), where 7 ' is the strength of the 
RPA perturbation. Note that at each time t used 
for solving the TDsBSE, one has to apply a differ¬ 
ent perturbation v'^ (r, t) at t = 0 , which is used to 
indicate the time for the RPA propagation. 


3. Propagate the orbitals using the adiabatic stochas- 



1 

2 


d^rd^r'lp^ {r,r',t)f 


xW^^^{r,r',0){r-r'). (31) 


In the above, the commuter vanishes for the k^x (i) term 
due to symmetry, but there is no a-priori reason why the 
k^x (t) term should vanish. However, as illustrated nu¬ 
merically in Fig.[^ the contribution of this non-vanishing 
term is rather small even on timescales much larger than 
the typical frequency in the system. 




















10 


^ S. Coe, W.-K. Woo, and V. B. Moungi Bawendi, Nature 
420, 800 (2002). 

^ N. Tessler, V. Medvedev, M. Kazes, S. H. Kan, and 
U. Banin, Science 295, 1506 (2002). 

^ I. Gur, N. A. Fromer, M. L. Geier, and A. P. Alivisatos, 
Science 310, 462 (2005). 

^ D. V. Talapin and C. B. Murray, Science 310, 86 (2005). 
® E. Runge and E. K. U. Gross, Phys. Rev. Lett. 52, 997 
(1984). 

® R. van Leeuwen, Inter. J. Moder. Phys. B 15, 1969 (2001). 
^ G. Onida, L. Reining, and A. Rubio, Rev. Mod. Phys. 
74, 601 (2002). 

® N. T. Maitra, K. Burke, H. Appel, and E. K. U. Gross, 
“Ten topical questions in time dependent density func¬ 
tional theory,” in Reviews in Modem Quantum Chemistry: 
A celebration of the contributions of R. G. Parr, Vol. II, 
edited by K. D. Sen (World-Scientific, Singapore, 2002) 

p. 1186. 

® M. Marques and E. Gross, Annu. Rev. Phys. Ghem. 55, 
427 (2004). 

K. Burke, J. Werschnik, and E. K. U. Gross, J. Ghem. 
Phys. 123, 062206 (2005). 

S. Botti, A. Schindlmayr, R. Del Sole, and L. Reining, 
Rep. Prog. Phys. 70, 357 (2007). 

D. Jacquemin, E. A. Perpete, I. Ciofini, and C. Adamo, 
Acc. Ghem. Res. 42, 326 (2009). 

M. E. Casida, J. Mol. Struct. 914, 3 (2009). 

C. Adamo and D. Jacquemin, Ghem. Soc. Rev. 42, 845 
(2013). 

R. E. Stratmann, G. E. Scuseria, and M. J. Frisch, J. 
Ghem. Phys. 109, 8218 (1998). 

K. Yabana and G. F. Bertsch, Phys. Rev. B 54, 4484 
(1996). 

G. F. Bertsch, J. I. Iwata, A. Rubio, and K. Yabana, 
Phys. Rev. B 62, 7998 (2000). 

R. Baer and D. Neuhauser, J. Ghem. Phys. 121, 9803 
(2004). 

R. Bauernschmitt and R. Ahlrichs, Ghem. Phys. Lett. 
256, 454 (1996). 

R. Bauernschmitt, R. Ahlrichs, F. H. Hennrich, and 
M. M. Kappes, J. Am. Ghem. Soc. 120, 5052 (1998). 

J. Fabian, Theor. Ghem. Acc. 106, 199 (2001). 

I. Vasiliev, S. Ogut, and J. R. Chelikowsky, Phys. Rev. 
B 65, 115416 (2002). 

Y. H. Shao, M. Head-Gordon, and A. 1. Krylov, J. Ghem. 
Phys. 118, 4807 (2003). 

M. G. Troparevsky, L. Kronik, and J. R. Ghelikowsky, J. 
Ghem. Phys. 119, 2284 (2003). 

N. T. Maitra, J. Ghem. Phys. 122, 234104 (2005). 

J. Andzelm, A. M. Rawlett, J. A. Orlicki, and J. F. Sny¬ 
der, J. Ghem. Theory Comput. 3, 870 (2007). 

N. Govind, M. Valiev, L. Jensen, and K. Kowalski, J. 
Phys. Ghem. A 113, 6041 (2009). 

M. J. G. Peach, C. R. Le Sueur, K. Ruud, M. Guil¬ 
laume, and D. J. Tozer, Phys. Chem. Chem. Phys. 11, 
4465 (2009). 

N. Kuritz, T. Stein, R. Baer, and L. Kronik, J. Chem. 
Theory Comput. 7, 2408 (2011). 

M. J. G. Peach, M. J. Williamson, and D. J. Tozer, J. 
Ghem. Theory Comput. 7, 3578 (2011). 


M. Srebro, N. Govind, W. A. de Jong, and J. Autschbach, 
J. Phys. Chem. A 115, 10930 (2011). 

R. M. Richard and J. M. Herbert, J. Chem. Theory Com¬ 
put. 7, 1296 (2011). 

A. Chantzis, A. D. Laurent, C. Adamo, and 

D. Jacquemin, J. Chem. Theory Comput. 9, 4517 (2013). 
R. Bauernschmitt, M. Haser, O. Treutler, and 

R. Ahlrichs, Chem. Phys. Lett. 264, 573 (1997). 

J. R. Chelikowsky, L. Kronik, and 1. Vasiliev, J. Phys. 
Condes. Matrer 15, R1517 (2003). 

J. Gavnholt, A. Rubio, T. Olsen, K. S. Thygesen, and 
J. Schiotz, Phys. Rev. B 79, 195405 (2009). 

S. Hirata and M. Head-Gordon, Chem. Phys. Lett. 302, 
375 (1999). 

S. Hirata, T. J. Lee, and M. Head-Gordon, J. Ghem. 
Phys. Ill, 8904 (1999). 

D. Jacquemin, E. A. Perpete, G. E. Scuseria, 1. Giofini, 
and C. Adamo, J. Chem. Theory Comput. 4, 123 (2008). 

T. Stein, L. Kronik, and R. Baer, J. Chem. Phys. 131, 
244119 (2009). 

T. Stein, L. Kronik, and R. Baer, J. Am. Chem. Soc. 
131, 2818 (2009). 

H. Phillips, S. Zheng, A. Hyla, R. Laine, T. Goodson, 

E. Geva, and B. D. Dunietz, J. Phys. Chem. A 116, 1137 

( 2011 ). 

Q. Ou, S. Fatehi, E. Alguire, Y. Shao, and J. E. Subotnik, 

J. Chem. Phys. 141, 024114 (2014) 

M. Parac and S. Grimme, Chem. Phys. 292, 11 (2003). 

S. Grimme and M. Parac, Chem. Phys. Chem. 4, 292 
(2003). 

A. Dreuw and M. Head-Gordon, J. Am. Ghem. Soc. 126, 
4007 (2004). 

N. T. Maitra, F. Zhang, R. J. Cave, and K. Burke, J. 
Chem. Phys. 120, 5932 (2004). 

W. Hieringer and A. Gorling, Chem. Phys. Lett. 419, 557 
(2006). 

B. G. Levine, G. Ko, J. Quenneville, and T. J. Martinez, 
Mol. Phys. 104, 1039 (2006). 

K. Lopata, R. Reslan, M. Kowaska, D. Neuhauser, 
N. Govind, and K. Kowalski, J. Chem. Theory Comput. 
7, 3686 (2011). 

T. Kowalczyk, S. R. Yost, and T. Van Voorhis, J. Chem. 
Phys. 134, 054128 (2011) 

C. M. Isborn, B. D. Mar, B. F. Curchod, 1. Tavernelli, 
and T. J. Martinez, J. Phys. Chem. B 117, 12189 (2013) 
S. Albrecht, L. Reining, R. Del Sole, and G. Onida, Phys. 
Rev. Lett. 80, 4510 (1998). 

M. Rohlfing and S. G. Louie, Phys. Rev. B 62, 4927 

( 2000 ). 

F. Sottile, M. Marsili, V. Olevano, and L. Reining, Phys. 
Rev. B 76, 161103 (2007). 

L. Ramos, J. Paier, G. Kresse, and F. Bechstedt, Phys. 
Rev. B 78, 195423 (2008). 

D. Rocca, Y. Ping, R. Gebauer, and G. Galli, Phys. Rev. 
B 85, 045116 (2012). 

L. Hedin, Phys. Rev. 139, A796 (1965). 

E. E. Salpeter and H. A. Bethe, Phys. Rev. 84, 1232 
(1951). 

M. S. Hybertsen and S. G. Louie, Phys. Rev. Lett. 55, 
1418 (1985). 





11 


M. S. Hybertsen and S. G. Louie, Phys. Rev. B 34, 5390 
(1986). 

L. Steinbeck, A. Rubio, L. Reining, M. Torrent, 1. White, 
and R. Godby, Gomput. Phys. Gommun. 125, 05 (1999). 

M. M. Rieger, L. Steinbeck, I. White, H. Rojas, and 

R. Godby, Gomput. Phys. Gommun. 117, 211 (1999). 

®'* P. Rinke, A. Qteish, J. Neugebauer, C. Freysoldt, and 

M. Scheffler, New J. Phys. 7, (2005). 

®® J. B. Neaton, M. S. Hybertsen, and S. G. Louie, Phys. 
Rev. Lett. 97, 216405 (2006). 

®® M. L. Tiago and J. R. Ghelikowsky, Phys. Rev. B 73, 
205334 (2006). 

®^ G. Friedrich and A. Schindlmayr, NIG Series 31, 335 
(2006). 

®® M. Gruning, A. Marini, and A. Rubio, J. Ghem. Phys. 
124, 154108 (2006). 

®® M. Shishkin and G. Kresse, Phys. Rev. B 75, 235102 
(2007). 

P. Huang and E. A. Garter, Annu. Rev. Phys. Ghem. 59, 
261 (2008). 

C. Rostgaard, K. W. Jacobsen, and K. S. Thygesen, Phys. 
Rev. B 81, 085103 (2010). 

I. Tamblyn, P. Darancet, S. Y. Quek, S. A. Bonev, and 

J. B. Neaton, Phys. Rev. B 84, 201402 (2011). 

P. Liao and E. A. Garter, Phys. Ghem. Ghem. Phys. 13, 
15189 (2011). 

S. Refaely-Abramson, R. Baer, and L. Kronik, Phys. Rev. 
B 84, 075144 (2011). 

N. Marom, F. Garuso, X. Ren, O. T. Hofmann, 

T. Korzdorfer, J. R. Ghelikowsky, A. Rubio, M. Scheffler, 
and P. Rinke, Phys. Rev. B 86, 245127 (2012). 

^® L. Y. Isseroff and E. A. Garter, Phys. Rev. B 85, 235142 

( 2012 ). 

S. Refaely-Abramson, S. Sharifzadeh, N. Govind, 
J. Autschbach, J. B. Neaton, R. Baer, and L. Kronik, 
Phys. Rev. Lett. 109, 226405 (2012). 

L. Kronik, T. Stein, S. Refaely-Abramson, and R. Baer, 
J. Ghem. Theory Gomput. 8, 1515 (2012). 

L. X. Benedict, E. L. Shirley, and R. B. Bohn, Phys. Rev. 
Lett. 80, 4514 (1998). 

®® L. X. Benedict, A. Puzder, A. J. Williamson, J. G. Gross- 
man, G. Galli, J. E. Klepeis, J.-Y. Raty, and O. Pankra¬ 
tov, Phys. Rev. B 68, 085310 (2003). 

C. D. Spataru, S. Ismail-Beigi, L. X. Benedict, and S. G. 
Louie, Phys. Rev. Lett. 92, 077402 (2004). 

N. Sai, M. L. Tiago, J. R. Ghelikowsky, and F. A. Re- 
boredo, Phys. Rev. B 77, 161306 (2008). 

F. Fuchs, G. Rodl, A. Schleife, and F. Bechstedt, Phys. 
Rev. B 78, 085103 (2008). 

M. Palummo, G. Hogan, F. Sottile, P. Bagala, and A. Ru¬ 
bio, J. Ghem. Phys. 131, 084102 (2009). 

L. Schimka, J. Harl, A. Stroppa, A. Gruneis, M. Mars- 
man, F. Mittendorfer, and G. Kresse, Nat. Mater. 9, 741 
( 2010 ). 

®® D. Rocca, D. Lu, and G. Galli, J. Ghem. Phys. 133, 
164109 (2010). 

X. Blase and C. Attaccalite, Appl. Phys. Lett. 99, 171909 

( 2011 ). 

C. Faber, 1. Duchemin, T. Deutsch, C. Attaccalite, V. 01- 
evano, and X. Blase, J. Mater. Sci. 47, 7472 (2012). 

G. Faber, P. Boulanger, G. Attaccalite, 1. Duchemin, and 
X. Blase, Philos. Trans. A Math. Phys. Eng. Sci. 372, 
20130271 (2014). 


B. Walker, A. M. Saitta, R. Gebauer, and S. Baroni, 
Phys. Rev. Lett. 96, 113001 (2006) 

D. Rocca, R. Gebauer, Y. Saad, and S. Baroni, J. Ghem. 
Phys. 128, 154105 (2008). 

H. E. Wilson, F. Gygi, and G. Galli, Phys. Rev. B 78, 
113303 (2008) 

S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Gian- 
nozzi. Rev. Mod. Phys. 73, 515 (2001). 

R. Baer, D. Neuhauser, and E. Rabani, Phys. Rev. Lett. 
Ill, 106402 (2013) 

D. Neuhauser, R. Baer, and E. Rabani, J. Ghem. Phys. 
141, 041102 (2014). 

D. Neuhauser, E. Rabani, and R. Baer, J. Ghem. Theory 
Gomput. 9, 24 (2013). 

Q. Ge, Y. Gao, R. Baer, E. Rabani, and D. Neuhauser, 
J. Phys. Ghem. Lett. 5, 185 (2013). 

D. Neuhauser, E. Rabani, and R. Baer, J. Phys. Ghem. 
Lett. 4, 1172 (2013). 

R. Baer and E. Rabani, Nano Lett. 12, 2123 (2012). 

D. Neuhauser, Y. Gao, G. Arntsen, G. Karshenas, E. Ra¬ 
bani, and R. Baer, Phys. Rev. Lett. 113, 076402 (2014). 
Y. Gao, D. Neuhauser, R. Baer, and E. Rabani, J. Ghem. 
Phys. 142, 034106 (2015). 

M. E. Casida, Recent advances in density functional meth¬ 
ods 1, 155 (1995). 

M. E. Gasida, “Time-dependent density functional re¬ 
sponse theory of molecular systems: theory, computa¬ 
tional methods, and functionals,” in Recent Developments 
and Applications in Density Functional Theory, edited by 
J. M. Seminario (Elsevier, Amsterdam, 1996) pp. 391-439. 
F. Furche, Phys. Rev. B 64, 195120 (2001). 

F. J. Dyson, Phys. Rev. 90, 994 (1953). 

J. Taylor, Phys. Rev. 95, 1313 (1954). 

S. Hirata and M. Head-Gordon, Ghem. Phys. Lett. 314, 
291 (1999). 

S. Hirata, M. Head-Gordon, and R. J. Bartlett, J. Ghem. 
Phys. Ill, 10774 (1999). 

R. Baer and D. Neuhauser, Phys. Rev. Lett. 94, 043002 
(2005). 

E. N. Brothers, A. F. Izmaylov, J. O. Normand, 
V. Barone, and G. E. Scuseria, J. Ghem. Phys. 129, 
011102 (2008). 

L. W. Wang and A. Zunger, J. Phys. Ghem. 98, 2158 
(1994). 

L. W. Wang and A. Zunger, Phys. Rev. B 51, 17398 
(1995). 

L. W. Wang and A. Zunger, Phys. Rev. B 53, 9579 (1996). 
H. Fu and A. Zunger, Phys. Rev. B 55, 1642 (1997). 

A. J. Williamson and A. Zunger, Phys. Rev. B 61, 1978 

( 2000 ). 

A. Franceschetti and A. Zunger, Phys. Rev. B 62, 2614 

( 2000 ). 

A. Zunger, Physica Status Solidi B-Basic Research 224, 
727 (2001). 

H. X. Fu and A. Zunger, Phys. Rev. B 56, 1496 (1997). 

E. Rabani, B. Hetenyi, B. J. Berne, and L. E. Brus, J. 
Ghem. Phys. 110, 5355 (1999). 

F. A. Reboredo, A. Franceschetti, and A. Zunger, Phys. 
Rev. B 61, 13073 (2000). 

A. Franceschetti and A. Zunger, Phys. Rev. B 62, R16287 

( 2000 ). 

H. Eshet, M. Griinwald, and E. Rabani, Nano Lett. 13, 
5880 (2013). 


90 

91 

92 

y3 

94 

yo 

96 

97 

98 

99 

100 

101 

102 

103 

104 

105 

106 

107 

108 

109 

110 

111 

112 

113 

114 

115 

116 

117 

118 

119 

120 

121 

122 



12 


123 

124 


L. W. Wang and A. Zunger, Phys. Rev. Lett. 73, 1039 R. Kosloff, J. Phys. Chem. 92, 2087 (1988). 

(1994). 

A. Zunger and L. W. Wang, Appl. Surf. Sci. 102, 350 
(1996). 



